% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%latest date 30 OCT 2005
clear;
data = load('c:\manoj\projects\eej\eej_par_all.txt');
load c:\manoj\projects\eej\OBSdH1;
load c:\manoj\projects\eej\mjdkpf107 mjd kp;
%-------- correlate the CHAMP derived EEJ strength and ground delta H
% with LT correction
%--------------------------------------------------------------------

bin5 = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
bin10 = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
bin20 = [-60 -40 -20 0 20 40 60];

%Control Center

obsdata = 6;  %ETT-HYB = 1,ETT-ABG = 2;TIR-HYB = 3;TIR-ABG = 4; HUA-FUQ 5; MBO-GUI = 6; AAE-QSB 7; GUA-CBI 8; YAP-BIK 9;PND-HYB 10;
satdata = 2; % %Pear Current Density = 1; av cur distan 2;
obslt = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat = 10;
x1data = bin10;
mindate = tmjd(2000,01,01,0,0,0,1); % starting date
maxdate = tmjd(2002,12,31,0,0,0,1)+1; % ending date
sqcor = 0; %SQ correction
ltcor = 1; %LT correction
robft = 0; %Robust outlier elimination
plotopt = 'k--';
plotcol = 'k';
plotchar = 'c';

%load c:\manoj\projects\eej\guacm4.mat index sq_est

switch obsdata,
    case 1,%ETT-HYB
        k = nanmean(squeeze(OBSdH(2,:,:)-OBSdH(3,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(2,:,:)-OBSdH(3,:,:))',[1,1096*24]);
        load c:\manoj\projects\eej\ettcm4.mat index sq_est;
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5; %Oservatory Longitude of ETT is given as 82.5 as the IST is referenced to this longitude
        st = 12.8; % The bins are still centered around ETT 
        en = 132.8;
        plotoptions = 'b*-';

    case 2,%ETT-ABG
        k = nanmean(squeeze(OBSdH(2,:,:)-OBSdH(1,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(2,:,:)-OBSdH(1,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((78-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 3,%TIR-HYB
        k = nanmean(squeeze(OBSdH(4,:,:)-OBSdH(3,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(4,:,:)-OBSdH(3,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((77.85-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78.5 (TIR) and 82.5 (Allahabad - IST)
        xtime = [1:24];
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'b*-';

    case 4,%TIR-ABG
        k = nanmean(squeeze(OBSdH(4,:,:)-OBSdH(1,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2; 
        data1(:,1) = mjd;
        load c:\manoj\projects\eej\tircm4.mat index sq_est;
        data1(:,2) = reshape(squeeze(OBSdH(4,:,:)-OBSdH(1,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((75-82.5)/15)/24; % This assumes 75 (between TIR and ABG)
        xtime = [1:24];
        obslon = 82.5;
        st = 12.8;
        en = 132.8;
        plotoptions = 'k.-';

    case 5,%HUA-FUQ
        k = nanmean(squeeze(OBSdH(5,:,:)-OBSdH(6,:,:)),1);
        k = circshift(k,[0,19]);
        xtime = 0.5:1:23.5;       
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(5,:,:)-OBSdH(6,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obslon = 285;
%        st = -142;
%        en = -8;
        st  = -139.5
        en = -10.5
        load c:\manoj\projects\eej\huacm4.mat index sq_est;
        plotoptions = 'k*-';
        
    case 6,%MBO-GUI
        k = nanmean(squeeze(OBSdH(7,:,:)-OBSdH(8,:,:)),1);
        k = circshift(k,[0,23]); % 23 ~ = 343/15
        %xtime = 0.5:1:23.5;       
        xtime = 1:24;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(7,:,:)-OBSdH(8,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obs_index = 2;
        obslon = 343;
        st = -81.5;
        en = 47.5;
        plotoptions = 'g*-';
        load c:\manoj\projects\eej\mbocm4.mat index sq_est;
    case 7,%AAE-QSB
        k = nanmean(squeeze(OBSdH(9,:,:)-OBSdH(10,:,:)),1);
        k = circshift(k,[0,2]); % 2 ~ = 37.5/15
        xtime = 1:24;       
        obs_index = 2;   
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(9,:,:)-OBSdH(10,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        obslon = 37.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'k*-';
        load c:\manoj\projects\eej\aaecm4.mat index sq_est;
    case 8,%GUA-CBI
        k = nanmean(squeeze(OBSdH(11,:,:)-OBSdH(12,:,:)),1);
        k = circshift(k,[0,9]); % 9 ~ = 144/15 (the remaining would be taken care by shifting time axis to 1:24 (instead 0.5:1:23.5)
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(11,:,:)-OBSdH(12,:,:))',[1,1096*24]);
        data1(isnan(data1(:,2)),2) = 9999; 
        xtime = 1:24;       
        obs_index = 2;  
        obslon = 143;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'r*-';
        load c:\manoj\projects\eej\guacm4.mat index sq_est;
    case 9,%YAP=BIK
        load c:\manoj\projects\eej\YAPMNBIK;
        %load c:\manoj\projects\eej\YAPMNBIKHR;
        k = nanmean(YAPMNBIK);
        k = circshift(k,[1,554]); %554 = (138.5/15)*60 minutes (UT shift in minutes)
        %k = circshift(k,[1,15]); %554 = (138.5/15)*60 minutes (UT shift in minutes)
        xtime = (linspace(0,24,1440));       
        %xtime = 1:24;
        obs_index = 2;
        data1(:,1)= yapdate;
        %data1(:,1) = yaphrdate1;
        data1(:,2) = yapmnbik;
        %data1(:,2) = yapmnbikhr;
        %obslon = 138.5;
        obslon = 138.5;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'c*-';
   case 10,%PND-HYB
        k = nanmean(squeeze(OBSdH(13,:,:)-OBSdH(3,:,:)),1);
        k=circshift(k,[0,5]); % by circshifting the data 5 hoours and starting the center at 1 (not 0.5 hours) we achive UT to LT
        obs_index = 2;
        data1(:,1) = mjd;
        data1(:,2) = reshape(squeeze(OBSdH(13,:,:)-OBSdH(3,:,:))',[1,1096*24]);
        data1(:,1) = data1(:,1) + ((80-82.5)/15)/24; %NOTE this correction is for the Solar and IST difference between 78 (ETT) and 82.5 (Allahabad - IST)
        xtime = 1:24;
        obslon = 80; %Oservatory Longitude of ETT is given as 82.5 as the IST is referenced to this longitude
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        plotoptions = 'b*-';
     otherwise,
        display('Error');
        break;
        return;
end;


%INFORMATION
%The EEJ parameters in the file
%1-mjd,2-longitude of eq-crossinfg, 3-zero of current position-south
%4 -north,5-half-width of current - south,6-north
%7minimum current south,8-north,9-pear current density 10-average
%difference 11 - diff to south minumum 12 - diff to north minum, 13-total
%east current ,14 total return current., 15 -f,16 - eq angle, 17-
%eq-direction, 18 - mgnteic equator latitude
switch satdata,
    case 1,
        eej_index = 9;
    case 2,
        eej_index = 10;
    otherwise,
        display('Error');
        break;
        return;
end;        

if sqcor == 1,
    data1(index,obs_index) = data1(index,obs_index)-(sq_est(:,1)-sq_est(:,2));
   % data1(index,obs_index) = (sq_est(:,1)-sq_est(:,2));
end;

[p,s] = polyfit(xtime,k,9);%In case of minute data change the degree to 20
%----------------------------

data3(:,1)= data(:,1);
data3(:,2) = interp1(data1(:,1),data1(:,2),data(:,1));

ndd = zeros([1,25]);
nloop = 1;
lto = [];
ltc = [];
xdata = [];
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;

for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat & data(:,1) < maxdate  &  data(:,1) > mindate; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
ltt1 = [];
ltt2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data3(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data3(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    ltt1(i) = lt1;
    ltt2(i) = lt2;
    ltt3(i) = dummy1;
    ltt4(i) = dummy2;
    obst(i) = data3(ind,1);
    chpt(i) = IndiaE(i,1);
    
         if lt2 >= obslt(1) & lt2 <= obslt(2) & kp(ind) <= 2.5 & abs(data3(ind,obs_index)) < 300 & data3(ind,1) < max(data(:,1)), %The last fliter take care the nodata (9999) points
            ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,lt1); %Getting the expexted EEJ strength from poly coefficients for champ LT
            f_ob(ndata) = polyval(p,lt2); %Getting the expexted EEJ strength from poly coefficients for obs LT
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data3(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
    obst(ind2) = [];
    chpt(ind2) = [];
end;

if length(ind1) > 1,
    obsE = data3(ind1,:);
    X = IndiaE(:,eej_index);
    if ltcor == 1,
        Y = obsE(:,obs_index).*(f_ch./f_ob)';
    elseif ltcor == 0,
        Y = obsE(:,obs_index)';
    end;
    
    if robft == 1,    
%    Robust outlier rejection  -- START
     if length(X) > 3,
     [rs,rstat] = robustfit(X,Y);
     residuals = abs(rstat.resid);
     LLL = residuals == max(residuals); %Just remove the most noisy point
     X(LLL) = [];
     Y(LLL) = [];
     end;
    % Robust outlier rejection  --- END
    end;

    [cof,err] = corrcoef(X,Y);  % correlation & Correction for observatory data 
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    fprintf('%4.3f ',cof(2,1));
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata-1;
    lto(nloop) = min(ltt2);
    ltc(nloop) = max(ltt2);
    nloop = nloop +1;

    %----------- PLOT SCATTER
%subplot(2,3,nloop-1);
% if x1data(nbin) == 0,
% 
% % % figure(1);
% % % hold on;
% plot(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',[plotcol 's'],'MarkerFaceColor',plotcol);
% set(gca,'FontSize',16);
% grid on
% hold on;
% %axis([0 0.5 0 200]);
% h=xlabel('A/m');
% % %h=xlabel(['Longitude Bin ' sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ']);
% set(h,'FontSize',16);
% ylabel('\DeltaH');
% [c,sfit] = polyfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',1);
% m = polyval(c,[min(IndiaE(:,eej_index)),max(IndiaE(:,eej_index))]);
% % %
% % %plot(obst,chpt,'r.');
% % % % xcor1(nloop-1) =dummy(2,1);
% % [aaa,rstat] = robustfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)');% with added Zero the fitted curve is made to pass through origin
% % dummy = corrcoef(rstat.resid,(obst-chpt)*24);
%  plot([min(IndiaE(:,eej_index)),max(IndiaE(:,eej_index))],m,plotopt,'LineWidth', 2);
%  text(0,70,sprintf('dH = %5.1f * I + %5.2f',c(1),c(2)),'FontSize',16,'Color',plotcol);
% % end;
%  text(0,80,sprintf(' CC %4.2f',cof(2,1)),'FontSize',16);
% % %text(0,110,sprintf('dH = %5.1f * I + %5.2f',c(1),c(2)),'FontSize',16,'Color',plotcol);
% % h=title(sprintf('%d ',x1data(nbin)));
% % set(h,'FontSize',16);
% % %end;
% end;


%-----------------END
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt X Y;
end;

x1data(isnan(x1data)) = [];
fprintf('\n');
figure(2);
hold on;
errorbar(x1data,ceff,cerr/2,plotopt,'LineWidth',2);
axis([-80 80 -.6 1.4]);
grid on;
hold on;